v20200413.0726
Os exemplos aqui apresentados estão disponíveis. Caso queira usá-los, crie um projeto, coloque o arquivo desta aula e os seguintes arquivos na pasta do mesmo:
Correlação e regressão são as duas análises baseadas na distribuição multivariada. Uma distribuição multivariada é descrita como uma distribuição de múltiplas variáveis. A correlação é descrita como a análise que nos permite saber a associação ou a ausência do relacionamento entre duas variáveis ‘x’ e ‘y’. Por outro lado, a análise de regressão prediz o valor da variável dependente com base no valor conhecido da variável independente, assumindo que existe uma relação matemática média entre duas ou mais variáveis.
A diferença entre correlação e regressão é uma das perguntas mais frequentes em entrevistas. Muitas pessoas confundem os dois conceitos.
| Comparação por: | Correlação | Regressão |
|---|---|---|
|
Medida estatística que determina a associação entre duas variáveis | Descreve como uma variável independente é numericamente relacionada com a variável dependente |
|
Representação linear da relação entre duas variáveis | Ajuste da melhor função para estimar uma variável dependente com base na independente |
|
Adirecional. Não importa definir qual variável é independente ou dependente | Direcional, da variável independente para a dependente |
|
Adimensional. Indica em qual medida as duas variáveis variam conjuntamente | Dimensional. Indica o impacto da mudança de uma unidade de medida da variável conhecida (independente) na variável estimada (dependente) |
|
Encontrar um valor numérico para expressar a relação entre duas variáveis | Estimar o valor de uma variável aleatória com base nos valores de uma variável fixada |
Uma correlação é uma relação entre duas variáveis. Cada variável representa um fenômeno específico; quando um dos fenômenos muda em uma direção específica, a segunda variável muda na direção do primeiro ou na direção oposta do primeiro. A mudança dos dois fenômenos na mesma direção no sentido do aumento no primeiro deslocamento por um aumento no segundo ou vice-versa na redução do primeiro acompanhado por uma diminuição no segundo relacionamento é positiva ou crescente (associação positiva). Ao contrário, quando o aumento no primeiro deslocamento é associado com uma redução no segundo ou vice-versa, uma redução no primeiro fenômeno seja acompanhada por um aumento no segundo, dizemos que a ligação é inversa ou decrescente (associação negativa).
Uma regressão é um método pelo qual o valor de uma variável pode ser estimado pelo valor da outra variável pela equação de regressão. Tem os seguintes tipos:
Em tabelas de contingência, como acontece com o teste do qui-quadrado e em análise de exames diagnósticos, precisamos verificar a associação entre duas variáveis categóricas (nominais e ou ordinais).
No entanto, quando precisamos verificar a associação entre variáveis quantitativas uma tabela de contingência não pode mais nos atender, mesmo usando apenas variáveis discretas. Imagine, por exemplo, verificar associação entre idade e sono, ainda que computássemos com números inteiros, respectivamente em anos e horas, e estudássemos somente idosos. Não é possível criar uma categoria para cada idade indo de 60 a 100 anos e cada hora de sono indo de 2 a 10: a tabela teria 41 linhas x 9 colunas e muitas células vazias ainda que a amostra fosse de tamanho razoável.
Por exemplo, a partir dos dados obtidos de
Yunus RM, Wazid SW, Hairi NN, Choo WY, Hairi FM, Sooryanarayana R, et al. (2017) Association between elder abuse and poor sleep: A cross-sectional study among rural older Malaysians. PLoS ONE 12(7): e0180222. https://doi.org/10.1371/journal.pone.0180222
o seguinte RScript mostra o diagrama de dispersão (scatterplot) entre idades e escores do (PSQI) e constrói uma tabela de contingência da idade (linhas de 60 a 100 anos) contra o escore do questionário PSQI (que avalia qualidade do sono, colunas de 0 a 21):
source("tabelaPSQI.R")Número de indivíduos (idade x PSQI):
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21
60 - - 7 20 15 8 6 6 3 5 5 - 1 1 - - - - - - - -
61 - 1 17 17 27 27 15 5 5 3 7 1 1 - - - - - - - - -
62 - 2 8 26 26 11 11 10 5 1 5 2 - 1 - 1 1 - - - - -
63 - 3 6 20 18 6 9 7 4 3 2 4 - 1 - 1 - - - - - -
64 - - 10 30 34 15 10 5 6 2 - 1 1 2 - - 1 - - - - -
65 - 3 6 24 13 10 5 3 5 6 1 3 3 1 - - - - - - - -
66 - 2 8 13 25 11 6 4 7 3 1 - 3 2 - - - - - - - -
67 - 4 3 16 20 15 11 5 2 1 2 1 - 1 - 1 - - - - - -
68 - 1 3 5 14 9 7 1 - 2 3 2 1 - 1 - - - - - - -
69 1 1 3 9 11 7 5 2 3 - 3 - - - 1 - - - - - - -
70 - 3 5 13 12 9 3 3 1 3 1 1 2 1 - - - - - - - -
71 - - 4 11 23 5 6 5 3 4 1 1 - - - - - - - - - -
72 - 2 11 22 10 14 6 5 3 4 - 1 - 2 - - - - - - - -
73 - 1 6 20 15 11 8 7 1 6 - 2 1 1 - - - - - - - -
74 - 1 5 23 12 12 7 4 1 3 1 1 1 1 - - - - - - - -
75 - 4 6 18 23 5 8 6 4 4 - - 1 - - - - - - - 1 -
76 - - 6 13 13 8 6 3 4 5 2 - 1 - - - - - - - - -
77 - - 7 20 10 7 3 5 6 2 1 1 - 1 - - - - 1 - - -
78 - - 5 5 13 11 8 2 1 1 1 - - - - - - - - - - -
79 - - 3 7 13 3 3 3 - 2 3 1 - 1 - - - - - - - -
80 - 1 2 4 5 4 7 2 3 2 1 3 1 - - - - - - - - -
81 - - 2 6 6 2 1 3 - 1 2 1 - - - - - - - - - -
82 - - 3 4 6 2 3 2 2 1 1 - 2 - - - - - - - - -
83 - - - 4 2 2 - 1 - 1 - - - - - - - - - - - -
84 - - 1 4 3 5 1 2 1 2 - - - - - 1 - - - - - -
85 - - - 3 1 - 2 1 1 - - - 1 - - - - - - - - -
86 - 2 - 1 - 2 1 1 - 1 - - 1 - - - - - - - - -
87 - - - 1 2 2 - - 1 - - - - - - - - - - - - -
88 - - 2 - 1 1 1 - - - - - - - - - - - - - - -
89 - - - - - - 1 - - - - - - - - - - - - - - -
90 - - - - - - - - - - - - - - - - - - - - - -
91 - - - - - - - - - - - - - - - - - - - - - -
92 - - - - - - - - - - - - - - - - - - - - - -
93 - - - - - - - - - - 1 - - - - - - - - - - -
94 - - - - - - - - - - - - - - - - - - - - - -
95 - - - - - - - - - - - - - - - - - - - - - -
96 - - - - - - - - - - - - - - - - - - - - - -
97 - - - - - - - - - - - - - - - - - - - - - -
98 - - - - - - - - - - - - - - - - - - - - - -
99 - - - - - - - - - - - - - - - - - - - - - -
100 - - - - 1 - - - - - - - - - - - - - - - - -
Ainda que este estudo conte com 1648 indivíduos estudados, a tabela tem com muitas células vazias. Além disto, idade (variável quantitativa) e o escore do PSQI (que também pode ser admitido como variável quantitativa) foram demovidos a variáveis qualitativas nominais, portanto perde-se informação.
As análises de correlação ou de regressão simples são formas mais adequadas para avaliar as associações entre variáveis quantitativas. Para tanto, como dizem os nomes das análises, primeiro verificaremos se há correlação (associação) entre as variáveis de interesse e, então, podemos estender para a regressão e prever o valor da variável dependente a partir da independente.
Apresentamos dados de um estudo realizado no Centro de Sáude Geraldo de Paula Souza, da Faculdade de Saúde Pública da USP, no qual gestantes foram avaliadas quanto à deficiência de ácido fólico e vitamina B12 no início do pré-natal e no 7º mês de gestação. Um grupo controle de mulheres saudáveis foi constituído para servir de referência para as medidas das vitaminas acima mencionadas. Durante o parto, foram obtidas amostras de sangue do recém-nascido para verificar a concentração de anticorpos contra a rubéola.
A planilha Gestantes.xlsx contém os seguintes campos (colunas):
| Variável | Descrição |
|---|---|
| Grupo | Gestantes ou Controle |
| NOME | Iniciais do nome das pacientes |
| IDADE | Idade das pacientes em anos |
| COR | Etnia das pacientes: br = branca; pd = mulata; am = amarela |
| HB | Hemoglobina (g/dL) no início do pré-natal |
| HT | Hematócrito (%) no início do pré-natal |
| HEM | Contagem de Hemácias (milhões/mm³) no início do pré-natal |
| LEUC | Contagem de Leucócitos (milhares/mm³) no início do pré-natal |
| RET | Reticulócitos (%) no início do pré-natal |
| PLAQUET | Plaquetas (milhares/mm³) no início do pré-natal |
| FOLICO | Concentração de Ácido Fólico (ng/mL) no início do pré-natal (valores abaixo de 3 ng/mL indicam deficiência de Ácido Fólico) |
| B12 | Concentração de Vitamina B12 (pg/mL) no início do pré-natal (valores abaixo de 400 pg/mL indicam deficiência de Vitamina B12) |
| FOL_7m | Concentração de Ácido Fólico (ng/mL) no sétimo mês de gestação |
| B12_7m | Concentração de Vitamina B12 (pg/mL) no sétimo mês de gestação |
| Hb_7m | Hemoglobina no sétimo mês de gestação |
| Ht_7m | Hematócrito no sétimo mês de gestação |
| Hm_7m | Contagem de Hemácias no sétimo mês de gestação |
| Leu_7m | Contagem de Leucócitos no sétimo mês de gestação |
| Ret_7m | Reticulócitos no sétimo mês de gestação |
| Plq_7m | Contagem de Plaquetas no sétimo mês de gestação |
| AC_Rub_MAE | Anticorpos IgG contra o vírus da Rubéola da gestante (UI/mL) |
| AC_Rub_RN | Anticorpos IgG contra o vírus da Rubéola do recém-nascido (UI/mL) |
| Primigesta | Primeira gestação? sim ou não |
| TABAGISMO | Gestante fuma? sim ou não |
Há diversas variáveis quantitativas nestes dados, obtidos de 40 gestantes e 25 controles. Usaremos apenas algumas das variáveis
para responder à pergunta:
Podemos estimar a concentração de hemoglobina e hemograma (contagem de eritrócitos e leucócitos) no sangue pela medida do hematócrito?
A coleta é mais simples; o processamento requer uma centrífuga e uma régua:
A coleta do sangue é mais trabalhosa:
O equipamento para a medida é mais sofisticado; os mais simples precisam do equipamento e de um computador:
Os mais sofisticados automatizam o processo:
Retomando a planilha Gestantes.xlsx, começamos por uma análise descritiva executando Gestantes_descritiva.R:
source("Gestantes_descritiva.R")
Hematocrito:
Total:
Min. 1st Qu. Median Mean 3rd Qu. Max.
34.00 36.00 38.00 38.58 41.00 45.00
por grupo:
$controle
Min. 1st Qu. Median Mean 3rd Qu. Max.
36.00 39.00 41.00 40.48 43.00 44.00
$Gestantes
Min. 1st Qu. Median Mean 3rd Qu. Max.
34.0 36.0 37.0 37.4 39.0 45.0
Hemoglobina:
Total:
Min. 1st Qu. Median Mean 3rd Qu. Max.
11.20 12.50 13.30 13.43 14.20 16.40
por grupo:
$controle
Min. 1st Qu. Median Mean 3rd Qu. Max.
12.30 13.40 14.20 14.05 14.60 15.80
$Gestantes
Min. 1st Qu. Median Mean 3rd Qu. Max.
11.20 12.40 12.75 13.04 13.60 16.40
Hemacias:
Total:
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.530 4.240 4.440 4.498 4.650 5.830
por grupo:
$controle
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.180 4.430 4.650 4.785 4.910 5.830
$Gestantes
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.530 4.112 4.280 4.318 4.490 5.380
Leucocitos:
Total:
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.600 6.400 7.600 7.902 9.000 15.400
por grupo:
$controle
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.600 5.100 6.300 6.528 7.000 13.800
$Gestantes
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.20 7.60 8.40 8.76 9.70 15.40
É possível que exista uma relação entre hematócrito (HT) e concentração de hemoglobina (HB) e entre hematócrito e número de hemácias (HEM). Entre hematócrito e o número de leucócitos é duvidoso.
Além disto, hematócrito parece mais bem associado com hemoglobina do que com o número de hemácias pois ambas aumentam quando o hematócrito aumenta (relação positiva ou crescente); aparentemente, uma reta parece ser capaz de explicar o relacionamento entre estas variáveis (relação linear).
Precisamos verificar se nossas impressões se sustentam e se podem ser quantificadas.
O primeiro passo é verificar se as variáveis são correlacionadas. Não há sentido fazer regressão entre elas, para prever o valor médio da variável dependente a partir de uma variável independente conhecida quando sequer estejam relacionadas. Mais do que isto, para uma regressão linear, o coeficiente de correlação precisa ser estatisticamente diferente de zero.
|
É comum encontrarmos nos softwares estatísticos 3 coeficientes de correlação: Coeficiente de correlação de Pearson Coeficiente de correlação de Spearman \(\tau\) de Kendall |
Também chamado de:
Avalia o grau de linearidade de duas variáveis quantitativas sob as suposições de:
O Teste da Hipótese Nula e o Intervalo de Confiança estimados por fórmula necessitam das suposições de:
O Intervalo de Confiança estimado por reamostragem (bootstrapping) não necessita da suposição de normalidade bivariada.
O valor absoluto do \(r\) é a intensidade da linearidade; seu quadrado é a proporção da variância compartilhada. Ambos são medidas de tamanho de efeito:
Observe: o \(r\) de Pearson serve apenas para para associações LINEARES.
Admita duas variáveis quantitativas fictícias, em pares (\(x\),\(y\)):
source("friendlycolor.R")
# valores
x <- 2:10
y <- c(1, 2.5, 2.4, 11, 5.8, 13.8, 14, 8.3, 13)
# para nao ter um ponto na coordenada do centroide
# (melhora o exemplo)
centro_x <- mean(x)
x[x==centro_x] <- x[x==centro_x]-0.5
cat("x:",x,"\n")x: 2 3 4 5 5.5 7 8 9 10
cat("y:",y,"\n")y: 1 2.5 2.4 11 5.8 13.8 14 8.3 13
# scatterplot
plot(x, y,
xlim = c(0,12), ylim = c(0,16),
pch=21, col="black", bg=friendlycolor(24), lty=2
)Aparentemente há uma relação crescente de \(y\) em função de \(x\).
Caso esta relação não existisse, deveríamos esperar que \(y\) não variasse com \(x\), e a relação seria neutra, representada por uma reta horizontal. Os valores observados de \(y\) flutuariam ao redor do valor médio de \(y\) (linha horizontal pontilhada):
# media de y, traca linha horizontal
centro_y <- mean(y)
lines (x, rep(centro_y,length(x)), col=friendlycolor(30),
lwd=3, lty=2)Vamos, então, admitir que \(y\) siga uma função linear de \(x\), na forma \[y = \beta_0 + \beta_1 x\] onde: \[\beta_0 ~ \text{... intercepto}\] \[\beta_1 ~ \text{... inclinação}\]
Podemos facilmente aplicar um modelo linear com a função lm() do R, obtendo o intercepto e a inclinação da reta que melhor se ajusta aos dados observados \((x,y)\). Como já há vários elementos neste gráfico, adicionaremos também uma legenda, ecoando a equação da reta que encontramos:
# modelo linear, traca a reta
modelo <- lm(y ~ x)
intercepto <- modelo$coefficients[1]
inclinacao <- modelo$coefficients[2]
y.estimado <- intercepto + inclinacao*x
lines (x, y.estimado, col=friendlycolor(8),
lwd=3, lty=1)
# legenda
legend ("topleft",
c("y observados","y médio",
paste("reta: y_medio = ",round(intercepto,3)," + ",
round(inclinacao,3),"x",sep="")
),
lwd=c(NA, 3, 3),
pch=c(21, NA, NA),
col=c("black",friendlycolor(30), friendlycolor(8)),
pt.bg=c(friendlycolor(24),friendlycolor(30), friendlycolor(8)),
cex=0.8, box.lwd=0)Neste exemplo, a reta foi ajustada pelo método de mínimos quadrados.
Infinitas retas passam pelo centróide \((\bar{x},\bar{y})\), as coordenadas dos valores médios dos pares \((x,y)\) observados. Construindo-se quadrados cujos lados sejam a distância entre os valores observados de \(y\) e os estimados, a reta que melhor se ajusta é aquela cuja somatória das áreas dos quadrados for mínima.
Rode demo_minimos_quadrados.R para uma demonstração:
source("demo_minimos_quadrados.R")Esta reta ajustada também é a que minimiza o valor de \(r\), medida do melhor ajuste encontrado. Em R, calculamos \(r\) e, simultaneamente, testamos a hipótese nula da nulidade de \(\rho\) (do qual \(r\) é o estimador) para concluirmos se, populacionalmente, a correlação existe:
\[H_0: \rho = 0\] \[H_1: \rho \ne 0\] \[\alpha=0.05 ~~\text{- default da função R}\]
com:
x: 2 3 4 5 5.5 7 8 9 10
y: 1 2.5 2.4 11 5.8 13.8 14 8.3 13
correlacao <- cor.test(x, y)
print(correlacao)
Pearson's product-moment correlation
data: x and y
t = 3.5904, df = 7, p-value = 0.008852
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3028074 0.9573293
sample estimates:
cor
0.8050357
Uma das primeiras confusões entre correlação e regressão vem de exemplos como este que acabamos de usar. Foi dito, acima, que correlação é adimensional e adirecional, servindo para medir a intensidade da associação entre duas variáveis; regressão, ao contrário, é dimensional e direcional, servindo para, a partir de uma VI conhecida, estimar o valor médio de uma VD.
O exemplo que acabamos de ver iniciou-se pela reta de regressão para demonstrar a correlação, dando a entender que a correlação, portanto, é dimensional e direcional. Não é.
Apesar de utilizarmos cor.test() com os valores brutos, o \(r\) de Pearson é um coeficiente padronizado: internamente é computada a intensidade da associação das duas variáveis padronizadas. Padronizar é computar o escore z, subtraindo-se a média e dividindo-se todos os valores pelo desvio-padrão: \[z_x = {{x_i - \bar{x}}\over{s_x}}\] e \[z_y = {{y_i - \bar{y}}\over{s_y}}\] Os respectivos escores \(z\) são adimensionais e a forma das respectivas distribuições não é afetada.
Podemos mostrar, usando o R como um laboratório, que a correlação não é afetada pela padronização, mas a regressão difere. Compare as saídas de Correlacao_e_Regressao_Padronizadas.R:
# Correlacao_e_Regressao_Padronizadas.R
# valores
x <- 2:10
y <- c(1, 2.5, 2.4, 11, 5.8, 13.8, 14, 8.3, 13)
# para nao ter um ponto na coordenada do centroide\
# (melhora o exemplo)
centro_x <- mean(x)
x[x==centro_x] <- x[x==centro_x]-0.5
cat("\nValores originais:\n")
cat("x:",x,"\n")
cat("y:",y,"\n")
cat("\nCorrelação com as variáveis originais:\n")
correlacao <- cor.test(x, y)
print(correlacao)
cat("\nRegressão com as variáveis originais:\n")
modelo <- lm(y ~ x)
print (modelo)
# padronizando x e y
centro_x <- mean(x)
desvpad_x <- sd(x)
x_padrao <- (x-centro_x)/desvpad_x
centro_y <- mean(y)
desvpad_y <- sd(y)
y_padrao <- (y-centro_y)/desvpad_y
cat("\nValores padronizados:\n")
cat("x padronizada:",round(x_padrao,3),"\n")
cat("y padronizada:",round(y_padrao,3),"\n")
cat("\nCorrelação com as variáveis padronizadas:\n")
correlacao_padrao <- cor.test(x_padrao, y_padrao)
print(correlacao_padrao)
cat("\nRegressão com as variáveis padronizadas:\n")
modelo_padrao <- lm(y_padrao ~ x_padrao)
print (modelo_padrao)Obtendo-se as saídas:
Valores originais:
x: 2 3 4 5 5.5 7 8 9 10
y: 1 2.5 2.4 11 5.8 13.8 14 8.3 13
Correlação com as variáveis originais:
Pearson's product-moment correlation
data: x and y
t = 3.5904, df = 7, p-value = 0.008852
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3028074 0.9573293
sample estimates:
cor
0.8050357
Regressão com as variáveis originais:
Call:
lm(formula = y ~ x)
Coefficients:
(Intercept) x
-1.152 1.536
Valores padronizados:
x padronizada: -1.438 -1.073 -0.709 -0.344 -0.162 0.385 0.749 1.114 1.478
y padronizada: -1.333 -1.047 -1.066 0.577 -0.416 1.112 1.151 0.062 0.959
Correlação com as variáveis padronizadas:
Pearson's product-moment correlation
data: x_padrao and y_padrao
t = 3.5904, df = 7, p-value = 0.008852
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3028074 0.9573293
sample estimates:
cor
0.8050357
Regressão com as variáveis padronizadas:
Call:
lm(formula = y_padrao ~ x_padrao)
Coefficients:
(Intercept) x_padrao
0.000 0.805
Retomando o exemplo do r de Pearson com as variáveis padronizadas torna o entendimento do \(r\) menos abstrato.
source("friendlycolor.R")
# valores
x <- 2:10
y <- c(1, 2.5, 2.4, 11, 5.8, 13.8, 14, 8.3, 13)
# para nao ter um ponto na coordenada do centroide
# (melhora o exemplo)
centro_x <- mean(x)
x[x==centro_x] <- x[x==centro_x]-0.5
# padronizando as variáveis x e y
centro_x <- mean(x)
desvpad_x <- sd(x)
z_x <- (x-centro_x)/desvpad_x
centro_y <- mean(y)
desvpad_y <- sd(y)
z_y <- (y-centro_y)/desvpad_y
cat("\nValores padronizados:\n")
cat("z_x:",round(z_x,3),"\n")
cat("z_y:",round(z_y,3),"\n")
cat("\nCorrelação com as variáveis padronizadas:\n")
z_correlacao <- cor.test(z_x, z_y)
print(z_correlacao)
# scatterplot
plot(z_x, z_y,
xlim = c(-2,2), ylim = c(-2,2),
pch=21, col="black", bg=friendlycolor(24)
)
# media de y, traca linha horizontal
z_centro_x <- mean(z_x)
z_centro_y <- mean(z_y)
lines (z_x, rep(z_centro_y,length(z_x)), col=friendlycolor(30),
lwd=3, lty=2)
# modelo linear, traca a reta
z_modelo <- lm(z_y ~ z_x)
z_intercepto <- z_modelo$coefficients[1]
z_inclinacao <- z_modelo$coefficients[2]
z_y.estimado <- z_intercepto + z_inclinacao*z_x
lines (z_x, z_y.estimado, col=friendlycolor(8),
lwd=3, lty=1)
# centroide
points(z_centro_x, z_centro_y, pch=21, col="black", bg="black")
# legenda
legend ("topleft",
c("y obs. pad.","y médio pad.",
paste("reta: y_medio = ",round(z_intercepto,3)," + ",
round(z_inclinacao,3),"x",sep="")
),
lwd=c(NA, 3, 3),
pch=c(21, NA, NA),
col=c("black",friendlycolor(30), friendlycolor(8)),
pt.bg=c(friendlycolor(24),friendlycolor(30), friendlycolor(8)),
cex=0.8, box.lwd=0)
Valores padronizados:
z_x: -1.438 -1.073 -0.709 -0.344 -0.162 0.385 0.749 1.114 1.478
z_y: -1.333 -1.047 -1.066 0.577 -0.416 1.112 1.151 0.062 0.959
Correlação com as variáveis padronizadas:
Pearson's product-moment correlation
data: z_x and z_y
t = 3.5904, df = 7, p-value = 0.008852
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.3028074 0.9573293
sample estimates:
cor
0.8050357
Com a padronização não é surpresa que o centróide seja \((0,0)\) e, consequentemente, o intercepto seja nulo.
| O script completo, ilustração o \(r\) de Pearson, está disponível em demo_r_de_Pearson.R. |
O \(r\) de Pearson é uma mera medida da intensidade de associação entre duas variáveis e a regressão em si, faz pouco sentido. No entanto, visualizá-la pode ajudar: é possível treinar-se para acertar o valor de \(r\) observando um gráfico desde que as variáveis sejam padronizadas antes; caso contrário, a distorção causada pela unidades de medida das variáveis \(x\) e \(y\) atrapalharão.
|
Uma implementação do jogo está no APEx. A partir do menu principal acesse: e escolha o jogo r in R. |
Varia de -1 (correlação negativa perfeita) a +1 (correlação positiva perfeita), passando pelo 0 (ausência de correlação).
Graficamente, quanto mais bem alinhados estiverem os pontos, mais \(|r|\) se aproxima de 1:
| \(r = -1\) | \(r = -0.8\) | \(r = -0.4\) | \(r = -0.2\) |
|---|---|---|---|
|
|
|
|
|
| \(r = 0\) | |||
|
|
|||
| \(r = 0.2\) | \(r = 0.4\) | \(r = 0.8\) | \(r = 1\) |
|
|
|
|
|
Porém, cuidado com a inclinação nula:
| \(r = -1\) | \(r = -1\) | \(r = -1\) |
|---|---|---|
|
|
|
|
| \(r = 0\) | ||
|
|
||
| \(r = 1\) | \(r = 1\) | \(r = 1\) |
|
|
|
|
Rodando Gestantes_rPearson.R revisitamos os gráficos de hemoglobina, hemácias e leucócitos em função do hematócrito (agora padronizados), verificando suas respectivas correlações:
# Gestantes_rPearson.R
library(readxl)
source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")
col_HB <- eiras.friendlycolor(30) # tijolo
pch_HB <- 22
col_HEM <- eiras.friendlycolor(28) # bordo
pch_HEM <- 23
col_LEUC <- eiras.friendlycolor(9) # azul
pch_LEUC <- 24
Gestantes <- read_excel("Gestantes.xlsx")
# HT x HB
eiras.cor.test.boot (Gestantes$HT, Gestantes$HB,
alpha=0.05, B=0,
z.score=TRUE,
xlab="Hematocrito (z)", ylab="Hemoglobina (z)",
bg=col_HB, col="black", pch=pch_HB)
# HT x HEM
eiras.cor.test.boot (Gestantes$HT, Gestantes$HEM,
alpha=0.05, B=0,
z.score=TRUE,
xlab="Hematocrito (z)", ylab="Hemoglobina (z)",
bg=col_HEM, col="black", pch=pch_HEM)
# HT x LEUC
eiras.cor.test.boot (Gestantes$HT, Gestantes$LEUC,
alpha=0.05, B=0,
z.score=TRUE,
xlab="LEUCatocrito (z)", ylab="LEUCoglobina (z)",
bg=col_LEUC, col="black", pch=pch_LEUC)Obtendo:
Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
Hematocrito (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.5757 -0.8883 -0.2009 0.0000 0.8302 2.2050
Hemoglobina (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.0196 -0.8419 -0.1171 0.0000 0.6983 2.6914
Pearson's product-moment correlation
data: x and y
t = 15.477, df = 63, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.8249474 0.9315439
sample estimates:
cor
0.8898127
Hematocrito (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.5757 -0.8883 -0.2009 0.0000 0.8302 2.2050
Hemoglobina (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.0312 -0.5407 -0.1208 0.0000 0.3201 2.7973
Pearson's product-moment correlation
data: x and y
t = 7.7299, df = 63, p-value = 1.059e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5468409 0.8046606
sample estimates:
cor
0.6976867
Hematocrito (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.5757 -0.8883 -0.2009 0.0000 0.8302 2.2050
Hemoglobina (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.9037 -0.6645 -0.1335 0.0000 0.4861 3.3186
Pearson's product-moment correlation
data: x and y
t = -1.6304, df = 63, p-value = 0.108
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4242872 0.0448923
sample estimates:
cor
-0.2012099
Vemos, portanto, que a concentração de hemoglobina e a contagem de hemácias têm correlação significante com o hematócrito para \(\alpha=0.05\), mas a contagem de leucócitos não tem.
Elaboramos Gestantes_rPearson_boot_z.R e funções acompanhantes para obter um método mais robusto, aqui exibindo a versão com regressão sobre as variáveis padronizadas:
# Gestantes_rPearson_boot_z.R
# (versao padronizada)
alfa <- 0.05
B <- 1e4
library(readxl)
library(dplyr)
source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")
col_HB <- eiras.friendlycolor(30) # tijolo
pch_HB <- 22
col_HEM <- eiras.friendlycolor(28) # bordo
pch_HEM <- 23
col_LEUC <- eiras.friendlycolor(9) # azul
pch_LEUC <- 24
Gestantes <- read_excel("Gestantes.xlsx")
# HT x HB
eiras.cor.test.boot (Gestantes$HT, Gestantes$HB,
z.score = TRUE,
alpha=alfa, B=B,
xlab="Hematocrito (z)", ylab="Hemoglobina (z)",
bg=col_HB, col="black", pch=pch_HB)
# HT x HEM
eiras.cor.test.boot (Gestantes$HT, Gestantes$HEM,
z.score = TRUE,
alpha=alfa, B=B,
xlab="Hematocrito (z)", ylab="Hemacias (z)",
bg=col_HEM, col="black", pch=pch_HEM)
# HT x LEUC
eiras.cor.test.boot (Gestantes$HT, Gestantes$LEUC,
z.score = TRUE,
alpha=alfa, B=B,
xlab="Hematocrito (z)", ylab="Leucocitos (z)",
bg=col_LEUC, col="black", pch=pch_LEUC)Produzindo, para \(B=1e4\) iterações:
Hematocrito (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.5757 -0.8883 -0.2009 0.0000 0.8302 2.2050
Hemoglobina (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.0196 -0.8419 -0.1171 0.0000 0.6983 2.6914
Pearson's r (bootstrapping with 10000 iterations):
r: 0.8913601 [0.8362414, 0.9304591]
t: 15.6078, df: 63
p: 2.516685e-23
Pearson's product-moment correlation
data: x and y
t = 15.477, df = 63, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.8249474 0.9315439
sample estimates:
cor
0.8898127
Hematocrito (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.5757 -0.8883 -0.2009 0.0000 0.8302 2.2050
Hemacias (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-2.0312 -0.5407 -0.1208 0.0000 0.3201 2.7973
Pearson's r (bootstrapping with 10000 iterations):
r: 0.7012755 [0.5357378, 0.8131753]
t: 7.807935, df: 63
p: 7.734528e-11
Pearson's product-moment correlation
data: x and y
t = 7.7299, df = 63, p-value = 1.059e-10
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.5468409 0.8046606
sample estimates:
cor
0.6976867
Hematocrito (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.5757 -0.8883 -0.2009 0.0000 0.8302 2.2050
Leucocitos (z)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-1.9037 -0.6645 -0.1335 0.0000 0.4861 3.3186
Pearson's r (bootstrapping with 10000 iterations):
r: -0.2045698 [-0.4053023, 0.0125995]
t: -1.658803, df: 63
p: 0.1020918
Pearson's product-moment correlation
data: x and y
t = -1.6304, df = 63, p-value = 0.108
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.4242872 0.0448923
sample estimates:
cor
-0.2012099
Compare os resultados convencionais com estes obtidos por bootstrapping. A conclusão é a mesma: somente a concentração de hemoglobina e a contagem de hemácias têm correlação com o hematócrito.
O coeficiente de correlação de Spearman, embora não tenha sido originalmente computado desta forma, coincide com o \(r\) de Pearson quando substituímos os valores numéricos por seus postos (rank). Embora exista a função rank() no R, este trabalho não é necessário porque cor.test() faz o ajuste, automaticamente, quando method=“spearman”.
O \(r\) de Spearman também serve para variáveis ordinais. Neste caso, a representação gráfica não tem sentido.
Da mesma forma que o \(r\) de Pearson, o \(r\) de Spearman varia de \(-1\) a \(+1\).
A informação provida por este coeficiente é sobre quanto o (de)crescimento dos valores é monotônico (i.e., VE cresce sempre ou decresce sempre quando a VI aumenta). Assim:
\[H_0:~\text{as duas variáveis não têm relação de (de)crescimento}\] \[H_1:~\text{as duas variáveis têm relação de (de)crescimento}\]
O exemplo hipotético que executamos acima, com o método de Spearman, mostra:
# demo_r_de_Spearman_boot.R
alfa <- 0.05
B <- 1e4
library(readxl)
library(dplyr)
source("eiras.jitter.R")
source("eiras.friendlycolor.R")
source("eiras.cor.test.boot.R")
# cores
col <- eiras.friendlycolor(31) # preto
bg <- eiras.friendlycolor(24) # amarelo
pch <- 21 # circulo
# valores
x <- 2:10
y <- c(1, 2.5, 2.4, 11, 5.8, 13.8, 14, 8.3, 13)
# para nao ter um ponto na coordenada do centroide
# (melhora o exemplo)
centro_x <- mean(x)
x[x==centro_x] <- x[x==centro_x]-0.5
cat("x:",x,"\n")
cat("y:",y,"\n")
eiras.cor.test.boot (x, y,
alpha=alfa, B=B,
method="spearman",
xlab="x", ylab="y",
col=col, bg=bg, pch=pch)x: 2 3 4 5 5.5 7 8 9 10
y: 1 2.5 2.4 11 5.8 13.8 14 8.3 13
x
Min. 1st Qu. Median Mean 3rd Qu. Max.
2.000 4.000 5.500 5.944 8.000 10.000
y
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.000 2.500 8.300 7.978 13.000 14.000
Spearman's r (bootstrapping with 10000 iterations):
r: 0.7614679 [0.1965812, 1]
p: 0.01712763
Spearman's rank correlation rho
data: x and y
S = 28, p-value = 0.02139
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.7666667
O coeficiente de Spearman é máximo (\(r=1\)) quando o crescimento é monotónico. Por exemplo, executando demo_r_de_Spearman_boot_monotonico.R:
source("demo_r_de_Spearman_boot_monotonico.R")x: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
y: 1.972 2.316 10.182 15.456 21.021 35.486 39.048 39.215 50.328 69.527 144.079 266.425 351.394 496.108 693.31 894.604 1178.866 1208.258 1422.695 1441.131
x
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 5.75 10.50 10.50 15.25 20.00
y
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.972 31.870 106.803 419.071 743.633 1441.131
Spearman's r (bootstrapping with 10000 iterations):
r: 1 [1, 1]
p: 0
Spearman's rank correlation rho
data: x and y
S = 0, p-value = 5.976e-06
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
1
Basta que algum valor não seja maior que o anterior, para \(r < 1\). Para isto, modifico um valor do exemplo anterior com:
y[14] <- y[13]-100 # criando valorObtendo:
source("demo_r_de_Spearman_boot_quasemonotonico.R")x: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
y: 1.972 2.316 10.182 15.456 21.021 35.486 39.048 39.215 50.328 69.527 144.079 266.425 351.394 251.394 693.31 894.604 1178.866 1208.258 1422.695 1441.131
x
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 5.75 10.50 10.50 15.25 20.00
y
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.972 31.870 106.803 406.835 743.633 1441.131
Spearman's r (bootstrapping with 10000 iterations):
r: 0.9984814 [0.9590909, 1]
p: 4.056654e-24
Spearman's rank correlation rho
data: x and y
S = 6, p-value = 6.135e-06
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
0.9954887
O valor será de \(r=-1\) para o decrescimento monotônico:
source("demo_r_de_Spearman_boot_monotonicoD.R")x: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
y: 1697.175 1693.511 1692.384 1681.655 1662.215 1640.678 1595.231 1531.388 1461.494 1380.687 1334.24 1312.594 1228.563 1177.282 1149.491 1094.791 1090.062 795.446 648.886 267.313
x
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 5.75 10.50 10.50 15.25 20.00
y
Min. 1st Qu. Median Mean 3rd Qu. Max.
267.3 1135.8 1357.5 1306.8 1646.1 1697.2
Spearman's r (bootstrapping with 10000 iterations):
r: -1 [-1, -1]
p: 0
Spearman's rank correlation rho
data: x and y
S = 2660, p-value = 5.976e-06
alternative hypothesis: true rho is not equal to 0
sample estimates:
rho
-1
Exemplificamos em Gestantes_matrizcorrelacoes.R.R a montagem de uma matriz de correlações e sugerimos a função GGally::ggcorr() que auxilia com informação visual. O código é:
# Gestantes_matrizcorrelacoes.R
library(readxl)
library(GGally)
Gestantes <- read_excel("Gestantes.xlsx")
mc <- data.frame(Gestantes$IDADE,
Gestantes$HT,
Gestantes$HB,
Gestantes$HEM,
Gestantes$LEUC,
Gestantes$FOLICO,
Gestantes$B12)
names(mc) <- c("Idade","HT","HB","HEM","LEUC","FOLICO","B12")
print(head(mc))
cat("\n...\n")
print(tail(mc, addrownums = FALSE, n=2L))
cat("\nMatriz de correlacoes:\n")
print(cor(mc)) # matriz de correlacoes
# grafico da matriz
print(GGally::ggcorr(mc,
nbreaks = 6,
label = TRUE,
label_size = 4,
color = "#888888"))Resultando em:
Loading required package: ggplot2
Registered S3 methods overwritten by 'ggplot2':
method from
[.quosures rlang
c.quosures rlang
print.quosures rlang
Registered S3 method overwritten by 'GGally':
method from
+.gg ggplot2
Attaching package: 'GGally'
The following object is masked from 'package:dplyr':
nasa
Idade HT HB HEM LEUC FOLICO B12
1 24 40 14.0 4.53 8.3 11.6 324
2 34 36 11.2 4.89 11.0 12.8 740
3 27 34 12.2 4.15 8.7 13.9 360
4 37 34 11.9 3.65 8.4 13.3 291
5 26 37 13.2 4.32 10.0 6.6 303
6 28 39 13.6 4.57 7.1 8.0 292
...
Idade HT HB HEM LEUC FOLICO B12
64 19 42 14.9 4.84 4.9 4.4 405
65 32 42 15.0 5.40 4.9 3.5 454
Matriz de correlacoes:
Idade HT HB HEM LEUC
Idade 1.00000000 0.06403429 0.03309290 -0.02254294 0.2344499
HT 0.06403429 1.00000000 0.88981271 0.69768672 -0.2012099
HB 0.03309290 0.88981271 1.00000000 0.65385018 -0.1876529
HEM -0.02254294 0.69768672 0.65385018 1.00000000 -0.1614452
LEUC 0.23444985 -0.20120990 -0.18765292 -0.16144522 1.0000000
FOLICO 0.06404258 0.21815022 0.04858327 0.23278578 -0.1062217
B12 0.06685935 0.15093130 0.10786434 0.25578900 -0.2032199
FOLICO B12
Idade 0.06404258 0.06685935
HT 0.21815022 0.15093130
HB 0.04858327 0.10786434
HEM 0.23278578 0.25578900
LEUC -0.10622167 -0.20321995
FOLICO 1.00000000 0.40394434
B12 0.40394434 1.00000000
O primeiro passo a fazer é a estatística descritiva, tanto numérica quanto gráfica, recomendável antes de qualquer outro procedimento. O segundo é verificar se as duas variáveis estão correlacionadas.
A análise de regressão é uma extensão da correlação. A regressão linear só pode ser usada se os dados forem adequadamente lineares.
A pergunta é direcional: “Quanto muda a VD em média se a VE mudar?”
O objetivo, aqui, é chegar em uma equação de regressão que relacione as duas variáveis envolvidas.
Seus dados parecem seguir uma reta?
Quarteto de Anscombe é o nome dado a quatro conjuntos de dados que têm estatísticas descritivas quase idênticas (como a média e a variância), mas que têm distribuições muito diferentes e aparências muito distintas quando exibidos graficamente. Cada conjunto de dados consiste de onze pontos (x,y). Eles foram construídos em 1973 pelo estatístico Francis Anscombe, com o objetivo de demonstrar tanto a importância de se visualizar os dados antes de analisá-los, quanto o efeito dos outliers e outras observações influentes nas propriedades estatísticas. Ele descreveu o artigo como tendo a finalidade de combater a impressão entre os estatísticos de que “cálculos numéricos são exatos, mas gráficos são aproximados/grosseiros.”
Os quatro conjuntos resultam em:
source("Anscombe.R")
x
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.0 6.5 9.0 9.0 11.5 14.0
y
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.260 6.315 7.580 7.501 8.570 10.840
Pearson's product-moment correlation
data: x and y
t = 4.2415, df = 9, p-value = 0.00217
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4243912 0.9506933
sample estimates:
cor
0.8164205
x
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.0 6.5 9.0 9.0 11.5 14.0
y
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.100 6.695 8.140 7.501 8.950 9.260
Pearson's product-moment correlation
data: x and y
t = 4.2386, df = 9, p-value = 0.002179
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4239389 0.9506402
sample estimates:
cor
0.8162365
x
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.0 6.5 9.0 9.0 11.5 14.0
y
Min. 1st Qu. Median Mean 3rd Qu. Max.
5.39 6.25 7.11 7.50 7.98 12.74
Pearson's product-moment correlation
data: x and y
t = 4.2394, df = 9, p-value = 0.002176
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4240623 0.9506547
sample estimates:
cor
0.8162867
x
Min. 1st Qu. Median Mean 3rd Qu. Max.
8 8 8 9 8 19
y
Min. 1st Qu. Median Mean 3rd Qu. Max.
5.250 6.170 7.040 7.501 8.190 12.500
Pearson's product-moment correlation
data: x and y
t = 4.243, df = 9, p-value = 0.002165
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.4246394 0.9507224
sample estimates:
cor
0.8165214
Executando-se demo_regressao.R:
x: 2 3 4 5 5.5 7 8 9 10
y: 1 2.5 2.4 11 5.8 13.8 14 8.3 13
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-4.3705 -1.4952 -0.9557 2.8653 4.4727
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.1517 2.7730 -0.415 0.69035
x 1.5358 0.4277 3.590 0.00885 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 3.319 on 7 degrees of freedom
Multiple R-squared: 0.6481, Adjusted R-squared: 0.5978
F-statistic: 12.89 on 1 and 7 DF, p-value: 0.008852
Neste gráfico:
Na regressão linear simples (RLS) são apenas estas duas e, por isso, acompanha a análise de correlação. Em regressão univariada múltipla (RLM) várias VEs contribuem para prever o valor da VD. Por isso são sempre análises direcionais.
A linha horizontal (pontilhada, em vermelho) corresponde à hipótese nula, de ausência de correlação, \[H_0: r = 0\] pois é a reta esperada para valores da VD que não variam em função da VE. Pela análise de correlação, rejeitando \(H_0\) de acordo com cor.test(), podemos prosseguir com a RLS.
A reta de regressão (sólida, em azul) foi ajustada aos valores \(y_i\) observados (círculos) com os valores de intercepto e inclinação e cruza com a horizontal em seu centróide, representado por um círculo preto no gráfico, nas coordenadas \((\bar{x_i},\bar{y_i})\). Os coeficientes foram obtidos em demo_regressao.R, nas linhas
modelo <- lm (y ~x)
...
res <- summary(modelo)O intercepto fica em res$coefficients[1]=-1.1516605 (\(\beta_0\)) e a inclinação em res$coefficients[2]=1.5357934 (\(\beta_1\)). A equação de uma reta é
\[\hat{y} = \beta_0 + \beta_1 x\]
portanto a linha sólida do gráfico é :
\(\hat{y_i}\) \(=\) -1.1516605 \(+\) 1.5357934 \(x_i\)
Note o circunflexo sobre o valor \(\hat{y}\), indicando que este é o valor estimado da VD, portanto médio, a partir de \(x\), assumida como VE com valor conhecido.
|
Formalmente, como os valores de y foram estimados, há um erro, \(\epsilon\), a ser considerado: \[y = \beta_0 + \beta_1 x + \epsilon\] Sinonímia:
|
Como a reta sólida é inclinada, diferencia-se da horizontal \(H_0\) e, portanto, com a equação encontrada, dada a VE, calculamos VD. No entanto, observando-se os valores da amostra (círculos amarelos), vemos que nem tudo foi capturado pela RLS.
Podemos dizer que o distanciamento da reta inclinada em relação à horizontal é a parte “explicada” pelo modelo. O comprimento das linhas verticais com pontilhado mais fino, em azul, representam a porcentagem da variância explicada pelo modelo linear (ligadas ao coeficiente de determinação). As linhas pretas em pontilhados e traços vão dos valores observados até a reta inclinada, representando aquilo o modelo de RLS não capturou, os resíduos (estimativas das diferencas entre o valor observado da VD e o valor estimado pela regressão em cada ponto observado).
O coeficiente de determinação é dado por \(R^2\). Volte à saída textual do modelo linear e encontre
Multiple R-squared: 0.6481
Portanto, neste exemplo, \(R^2 = 64.81\%\) da variância amostral da VD foi explicada pela RLS. Os restantes \(35.19\%\) da variancia amostral da VD é explicada pelo termo de erro, i.e., por outras variáveis ou por outro formato de função (e.g., outro tipo de regressão não linear).
|
Em termos de analise de significancia pratica, os modelos lineares gerais, incluindo a analise de regressão, querem saber qual proporção da variância da VD é explicada pelos efeitos incluidos no modelo. Na RLS, há dois efeitos: a covariável VE e o termo de erro (tudo que não é VE). A estimativa amostral de \(R^2\) informa sobre a amostra. Para inferir em relação à população, usamos os intervalos de confiança do \(R^2\). |
Uma ligação entre a RLS e a correlação é que o valor do coeficiente de deteminação é numericamente coincidente com o quadrado do \(r\) de Pearson:
\[R^2 = r^2\]
Repare, em todos os gráficos anteriores, que as retas de regressão e respectivos intervalos de confiança não vão além dos valores de \(x_i\) observados. O motivo é que não sabemos se o fenômeno manterá o comportamento linear além deste ponto.
O terceiro conjunto de dados do Quarteto de Anscombe mostra este efeito, com um único ponto visivelmente deslocando a inclinação da reta de regressão. Dois exemplos, modificações de demo_regressao.R nas quais o ponto \((3, 2.4)\) foi deslocado para \((1.5, 13)\) e \((11,1)\):
source("demo_regressao_3.R")source("demo_regressao_3b.R")Para os métodos clássicos, a RLS requer:
Para os métodos mais robustos:
lm_robust() nao é um metodo por bootstrapping: dispensa apenas a violacao da suposicao de heterocedasticidade.
mcr::mcreg(), método por bootstrapping, só não lida com a independência entre dos pares.
Retomamos o exemplo para observar as premissas.
Já observamos nos gráficos de dispersão do hematócrito contra hemoglobina, hemácias e leucócitos. Parecem seguir uma reta uniformemente.
Em vez de testar rigorosamente a normalidade, vejamos seus density plots (linhas sólidas) e uma distribuição normal que tenha a mesma média e desvio-padrão dos dados (linhas pontilhadas).
Para o hematócrito:
library(readxl)
Gestantes <- read_excel("Gestantes.xlsx")
densidade <- density(Gestantes$HT)
media <- mean(Gestantes$HT)
desvpad <- sd(Gestantes$HT)
x <- seq(media-3*desvpad,
media+3*desvpad, length.out=100)
y <- dnorm(x, mean=media, sd=desvpad)
plot(densidade,
xlim=c(min(x),max(x)),
ylim=c(0, max(y,densidade$y,na.rm=TRUE)),
xlab="Hematocrito (%)",
ylab="densidade",
type="l")
lines(x,y,lty=2)concentração de hemoglobina:
library(readxl)
Gestantes <- read_excel("Gestantes.xlsx")
densidade <- density(Gestantes$HB)
media <- mean(Gestantes$HB)
desvpad <- sd(Gestantes$HB)
x <- seq(media-3*desvpad,
media+3*desvpad, length.out=100)
y <- dnorm(x, mean=media, sd=desvpad)
plot(densidade,
xlim=c(min(x),max(x)),
ylim=c(0, max(y,densidade$y,na.rm=TRUE)),
xlab="Hemoglobina (mg/dl)",
ylab="densidade",
type="l")
lines(x,y,lty=2)contagem de hemácias:
library(readxl)
Gestantes <- read_excel("Gestantes.xlsx")
densidade <- density(Gestantes$HEM)
media <- mean(Gestantes$HEM)
desvpad <- sd(Gestantes$HEM)
x <- seq(media-3*desvpad,
media+3*desvpad, length.out=100)
y <- dnorm(x, mean=media, sd=desvpad)
plot(densidade,
xlim=c(min(x),max(x)),
ylim=c(0, max(y,densidade$y,na.rm=TRUE)),
xlab="Hemácias (milhões/mm³)",
ylab="densidade",
type="l")
lines(x,y,lty=2)contagem de leucócitos:
library(readxl)
Gestantes <- read_excel("Gestantes.xlsx")
densidade <- density(Gestantes$LEUC)
media <- mean(Gestantes$LEUC)
desvpad <- sd(Gestantes$LEUC)
x <- seq(media-3*desvpad,
media+3*desvpad, length.out=100)
y <- dnorm(x, mean=media, sd=desvpad)
plot(densidade,
xlim=c(min(x),max(x)),
ylim=c(0, max(y,densidade$y,na.rm=TRUE)),
xlab="Leucócitos (milhares/mm³)",
ylab="densidade",
type="l")
lines(x,y,lty=2)Vamos assumir que as quatro variáveis têm distribuição amostral razoavelmente normal. Quanto à binormalidade (algo como um chapéu mexicano em 3 dimensões), não há teste disponível.
Montamos e generalizamos uma função em eiras.cor.test.boot.R que pode ser chamada para correlação e regressão, usando o método clássico ou bootstrapping, bastando alterar seus parâmetros. Foi utilizada a partir de Gestantes_rPearson.R (usando a função convencional do R) e Gestantes_rPearson_boot_z.R (usando bootstrapping), ambas com os dados padronizados e no contexto da correlação.
Para a RLS, o interesse em estabelecer uma equação para a previsão da VD necessita a utilização das unidades de medida. A mesma função, agora, será configurada para produzir os dados da regressão.
Como este exemplo aparenta ter linearidade, distribuição normal e homocedasticidade, utilizaremos os métodos clássicos, a partir de Gestantes_RLS.R:
|
Consulte o código:
|
source("Gestantes_RLS.R")
Hematocrito (%)
Min. 1st Qu. Median Mean 3rd Qu. Max.
34.00 36.00 38.00 38.58 41.00 45.00
Hemoglobina (mg/dl)
Min. 1st Qu. Median Mean 3rd Qu. Max.
11.20 12.50 13.30 13.43 14.20 16.40
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-1.39432 -0.26945 0.03055 0.31839 0.90568
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.40429 0.84391 0.479 0.634
x 0.33757 0.02181 15.477 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5077 on 63 degrees of freedom
Multiple R-squared: 0.7918, Adjusted R-squared: 0.7885
F-statistic: 239.5 on 1 and 63 DF, p-value: < 2.2e-16
Hematocrito (%)
Min. 1st Qu. Median Mean 3rd Qu. Max.
34.00 36.00 38.00 38.58 41.00 45.00
Hemacias (milhoes/mm³)
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.530 4.240 4.440 4.498 4.650 5.830
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-0.67232 -0.19498 -0.00498 0.14969 1.03346
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.09034 0.57174 0.158 0.875
x 0.11422 0.01478 7.730 1.06e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.3439 on 63 degrees of freedom
Multiple R-squared: 0.4868, Adjusted R-squared: 0.4786
F-statistic: 59.75 on 1 and 63 DF, p-value: 1.059e-10
Hematocrito (%)
Min. 1st Qu. Median Mean 3rd Qu. Max.
34.00 36.00 38.00 38.58 41.00 45.00
Leucocitos (milhares/mm³)
Min. 1st Qu. Median Mean 3rd Qu. Max.
3.600 6.400 7.600 7.902 9.000 15.400
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-4.5491 -1.4366 -0.2179 0.9759 7.0946
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 13.93071 3.70831 3.757 0.000379 ***
x -0.15626 0.09584 -1.630 0.108007
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 2.231 on 63 degrees of freedom
Multiple R-squared: 0.04049, Adjusted R-squared: 0.02526
F-statistic: 2.658 on 1 and 63 DF, p-value: 0.108
Observamos nas saídas as retas de regressão, a equação da reta, o coeficiente de determinação (\(R^2\)), e a banda de confiança de 95% (linhas tracejadas), onde esperamos que a reta de regressao populacional esteja contida.
O valor-\(p\) referente ao teste da hipótese nula:
\[H_0: \beta_1 = 0\] \[H_1: \beta_1 \ne 0\] Portanto, para \(\alpha=0.05\) só há modelo estatisticamente válido para estimar a concentração de hemoglobina e a contagem de hemácias a partir do hematócrito:
A impressão visual de que a predição para hemoglobina devia ser melhor que a para hemácias tem apoio na análise de regressão linear simples.
Ambos os tamanhos de efeito são grandes:
Ellis PD (2010) The Essential Guide to Effect Sizes, 1st ed., Cambridge University Press.
Este exemplo não atende às premissas da RLS clássica usada no exemplo anterior.
Podemos usar o tamanho do corpo para prever tamanho do cérebro em mamíferos a partir de alguns animais conhecidos?
Os dados estão em CorpoCerebro.xlsx. Rodamos a RLS clássica:
source("CorpoCerebro_RLS.R")
Corpo (kg)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 39.95 100.00 942.34 364.00 6654.00
Cerebro (g)
Min. 1st Qu. Median Mean 3rd Qu. Max.
5.5 168.0 419.0 1273.8 987.5 5712.0
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-692.8 -348.6 -240.3 -105.2 1887.7
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 427.2403 258.7882 1.651 0.133
x 0.8983 0.1200 7.484 3.76e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 772 on 9 degrees of freedom
Multiple R-squared: 0.8616, Adjusted R-squared: 0.8462
F-statistic: 56.01 on 1 and 9 DF, p-value: 3.755e-05
A saída textual parece mostrar tudo bem com a predição: \(R^2=86.16\%\) e \(p=3.76e-05\).
Lembre-se do quarteto de Anscombe!
Olhe os gráficos com atenção.
Os dados estão muito concentrados no lado inferior-esquerdo do gráfico (difícil aferir linearidade) e os dois elefantes afastados dos demais; os valores parecem tornar-se mais dispersos com os corpos maiores (sinal de heterocedasticidade). Além disto, as distribuições dos pesos dos corpos e dos cérebros não são aproximadamente normais:
Desenvolvemos uma pequena função capaz de testar as Transformações Potência de Tukey, indicando qual minimiza a assimetria (usando um coeficiente de assimetria também desenvolvido por Pearson). Está em eiras.tukey.try.R.
O peso do corpo:
library(readxl)
CorpoCerebro <- read_excel("CorpoCerebro.xlsx")
source("eiras.tukey.try.R")
potencia <- eiras.tukey.try(CorpoCerebro$Corpo, xlab="Corpo (kg)", show="all")
------------------
Tukey power = -3:
------------------
media: -0.09120363
desvio-padrao: 0.3014152
mediana: -1e-06
Coef. Assimetria de Pearson: -0.9077442
------------------
Tukey power = -2.5:
------------------
media: -0.09169442
desvio-padrao: 0.3012611
mediana: -1e-05
Coef. Assimetria de Pearson: -0.9130064
------------------
Tukey power = -2:
------------------
media: -0.09306477
desvio-padrao: 0.3008651
mediana: -1e-04
Coef. Assimetria de Pearson: -0.9269747
------------------
Tukey power = -1.5:
------------------
media: -0.0972522
desvio-padrao: 0.299867
mediana: -0.001
Coef. Assimetria de Pearson: -0.962949
------------------
Tukey power = -1:
------------------
media: -0.1128261
desvio-padrao: 0.2972685
mediana: -0.01
Coef. Assimetria de Pearson: -1.037709
------------------
Tukey power = -0.5:
------------------
media: -0.1961283
desvio-padrao: 0.2859996
mediana: -0.1
Coef. Assimetria de Pearson: -1.00834
------------------
Tukey power = 0:
------------------
media: 4.671906
desvio-padrao: 2.498658
mediana: 4.60517
Coef. Assimetria de Pearson: 0.08012583
------------------
Tukey power = 0.5:
------------------
media: 19.71794
desvio-padrao: 24.67576
mediana: 10
Coef. Assimetria de Pearson: 1.181476
------------------
Tukey power = 1:
------------------
media: 942.3364
desvio-padrao: 2033.755
mediana: 100
Coef. Assimetria de Pearson: 1.242534
------------------
Tukey power = 1.5:
------------------
media: 62798.11
desvio-padrao: 163680.2
mediana: 1000
Coef. Assimetria de Pearson: 1.132662
------------------
Tukey power = 2:
------------------
media: 4648142
desvio-padrao: 13284737
mediana: 10000
Coef. Assimetria de Pearson: 1.0474
------------------
Tukey power = 2.5:
------------------
media: 358772808
desvio-padrao: 1083301681
mediana: 1e+05
Coef. Assimetria de Pearson: 0.9932768
------------------
Tukey power = 3:
------------------
media: 28299251584
desvio-padrao: 88464183174
mediana: 1e+06
Coef. Assimetria de Pearson: 0.9596511
------------------
Melhor ajuste
------------------
logarithm transformation
Coef. Assimetria de Pearson: 0.08012583
deve ser transformado por logarítmo.
O peso do cérebro:
library(readxl)
CorpoCerebro <- read_excel("CorpoCerebro.xlsx")
source("eiras.tukey.try.R")
potencia <- eiras.tukey.try(CorpoCerebro$Cerebro, xlab="Cerebro (g)", show="y")
------------------
Melhor ajuste
------------------
logarithm transformation
Coef. Assimetria de Pearson: 0.1201111
também por logaritmo.
Então executamos CorpoCerebro_RLS_log.R, que cria duas colunas adicionais com a transformação sugerida:
CorpoCerebro$Corpo_log <- log(CorpoCerebro$Corpo)
CorpoCerebro$Cerebro_log <- log(CorpoCerebro$Cerebro)source("CorpoCerebro_RLS_log.R")
log[Corpo (kg)]
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 3.638 4.605 4.672 5.794 8.803
log[Cerebro (g)]
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.705 5.122 6.038 5.962 6.835 8.650
Call:
lm(formula = y ~ x)
Residuals:
Min 1Q Median 3Q Max
-1.0530 -0.4859 -0.2907 0.4574 1.5973
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.7577 0.5767 4.782 0.000999 ***
x 0.6858 0.1100 6.236 0.000152 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.8689 on 9 degrees of freedom
Multiple R-squared: 0.8121, Adjusted R-squared: 0.7912
F-statistic: 38.89 on 1 and 9 DF, p-value: 0.0001521
Os dados estão muito mais bem arranjados no gráfico. Para \(\alpha=0.05\) o modelo é estatisticamente válido (\(p=0.000152\)). O coeficiente de determinação é \(R^2=81.21\%\) (grande). Precisamos, no entanto, como usamos logarítmos, ajustar a equação:
\[ln(\hat{\text{Cérebro}}) = 2.7577 + 0.6858 \cdot ln(\text{Corpo})\] ou \[\hat{\text{Cérebro}} = e^{2.7577 + 0.6858 \cdot ln(\text{Corpo)}}\]
O lm_robust() só é robusto à heterocedasticidade. Então também pode ser executado com as transformações potência de Tukey:
source("CorpoCerebro_RLS_robust.R")Warning: Assigning non-quosure objects to quosure lists is deprecated as of rlang 0.3.0.
Please coerce to a bare list beforehand with `as.list()`
This warning is displayed once per session.
log[Corpo (kg)]
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.000 3.638 4.605 4.672 5.794 8.803
log[Cerebro (g)]
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.705 5.122 6.038 5.962 6.835 8.650
Call:
estimatr::lm_robust(formula = y ~ x)
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) 2.7577 0.7730 3.567 0.0060490 1.0090 4.5065 9
x 0.6858 0.1233 5.562 0.0003509 0.4069 0.9647 9
Multiple R-squared: 0.8121 , Adjusted R-squared: 0.7912
F-statistic: 30.94 on 1 and 9 DF, p-value: 0.0003509
desvio-padrao de log[Cerebro (g)] = 1.901556
Erro-padrao da estimativa de log[Cerebro (g)] = 0.8689489
R multiplo = 0.9011453
R^2 = coefic. de determinacao = 0.8120628
Eta^2 de log[Corpo (kg)] = 0.8120628
log[Cerebro (g)]_media = 2.757748 + 0.6857994 * log[Corpo (kg)] [0,8.802973]
log[Cerebro (g)]_media = 5.961738 + 0.6857994 * log[Corpo (kg)]_centrada [-4.671906,4.131068]
Não traz a banda de confiança, mas informa algumas estatísticas adicionais.
A função mcr::mcreg() é a mais robusta, por bootstrapping, e está exemplificada em RLS_robusta_Adm2008_Masc.R, prevendo a massa corpórea a partir da estatura:
source("RLS_robusta_Adm2008_Masc.R")Loading required package: car
Loading required package: carData
Attaching package: 'carData'
The following object is masked _by_ '.GlobalEnv':
Anscombe
Attaching package: 'car'
The following object is masked from 'package:dplyr':
recode
Loading required package: sandwich
Loading required package: tcltk
------------------------------------------
Reference method: Estatura (cm)
Test method: Massa Corporal Total (kg)
Number of data points: 51
------------------------------------------
The confidence intervals are calculated with
bootstrap ( quantile ) method.
Confidence level: 95%
------------------------------------------
LINEAR REGRESSION FIT:
EST SE LCI UCI
Intercept -96.1987808 NA -135.0318854 -55.584608
Slope 0.9677716 NA 0.7399102 1.190705
------------------------------------------
BOOTSTRAP SUMMARY
global.est bootstrap.mean bias bootstrap.se
Intercept -96.19878 -95.74980 0.44898 20.12221
Slope 0.96777 0.96535 -0.00242 0.11418
Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123
------------------------------------------
Reference method: Estatura (cm)
Test method: Massa Corporal Total (kg)
Number of data points: 51
------------------------------------------
The confidence intervals are calculated with
bootstrap ( quantile ) method.
Confidence level: 95%
------------------------------------------
LINEAR REGRESSION FIT:
EST SE LCI UCI
Intercept -96.1987808 NA -135.0318854 -55.584608
Slope 0.9677716 NA 0.7399102 1.190705
------------------------------------------
BOOTSTRAP SUMMARY
global.est bootstrap.mean bias bootstrap.se
Intercept -96.19878 -95.74980 0.44898 20.12221
Slope 0.96777 0.96535 -0.00242 0.11418
Bootstrap results generated with fixed RNG settings.
RNG kind: Mersenne-Twister
RNG seed: 123
mean sd 0% 25% 50% 75% 100% n
Estatura 176.35294 8.076691 156 172 176 182 193 51
MCT 74.47059 12.098517 48 66 73 82 105 51
NULL
Call:
estimatr::lm_robust(formula = MCT ~ Estatura, data = Dados)
Standard error type: HC2
Coefficients:
Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
(Intercept) -96.1988 19.6096 -4.906 1.068e-05 -135.6056 -56.792 49
Estatura 0.9678 0.1113 8.698 1.669e-11 0.7442 1.191 49
Multiple R-squared: 0.4174 , Adjusted R-squared: 0.4055
F-statistic: 75.66 on 1 and 49 DF, p-value: 1.669e-11
desvio-padrao de MCT = 8.076691
Erro-padrao da estimativa de MCT = 6.227405
R multiplo = 0.646062
R^2 = coefic. de determinacao = 0.4173961
Eta^2 de Estatura = 0.4173961
MCT.Media = -96.19878 + 0.9677716 * Estatura [ 156 , 193 ]
MCT.Media = 74.47059 + 0.9677716 * Estatura.Centrada [ -20.35294 , 16.64706 ]